Smooth k~temperature models instead of Sharpe-Schoolfield?

Load libraries

Code

pkgs <- c("tidyverse", "tidylog", "sdmTMB", "RColorBrewer", "viridis")

## minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.1.2     ✔ readr     2.1.4
#> ✔ forcats   1.0.0     ✔ stringr   1.5.0
#> ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
#> ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
#> ✔ purrr     1.0.1     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#> 
#> Attaching package: 'tidylog'
#> 
#> 
#> The following objects are masked from 'package:dplyr':
#> 
#>     add_count, add_tally, anti_join, count, distinct, distinct_all,
#>     distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
#>     full_join, group_by, group_by_all, group_by_at, group_by_if,
#>     inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
#>     relocate, rename, rename_all, rename_at, rename_if, rename_with,
#>     right_join, sample_frac, sample_n, select, select_all, select_at,
#>     select_if, semi_join, slice, slice_head, slice_max, slice_min,
#>     slice_sample, slice_tail, summarise, summarise_all, summarise_at,
#>     summarise_if, summarize, summarize_all, summarize_at, summarize_if,
#>     tally, top_frac, top_n, transmute, transmute_all, transmute_at,
#>     transmute_if, ungroup
#> 
#> 
#> The following objects are masked from 'package:tidyr':
#> 
#>     drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
#>     spread, uncount
#> 
#> 
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> 
#> 
#> Loading required package: viridisLite

# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Check the temperature model script! This is the order based on mean temperature, which makes sense for the sharpe schoolfield plot, and therefore we might keep it across plots
order <- data.frame(area = c("SI_HA", "BT", "TH", "SI_EK", "FM", "JM", "MU", "FB", "VN", "HO", "BS", "RA")) 

nareas <- length(unique(order))
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)

# Load functions
home <- here::here()
# fxn <- list.files(paste0(home, "/R/functions"))
# invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Read and trim data

Code

VBG <- read_csv(paste0(home, "/output/vbg.csv"))
#> Rows: 382 Columns: 10
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): area, type
#> dbl (8): cohort, k_mean, k_median, linf_median, k_lwr, k_upr, temp, n
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Fit smooth models of k vs temperature

Code

preds <- list()

for(a in unique(VBG$area)) {
  
  d <- filter(VBG, area == a)

  
  m <- sdmTMB(k_mean ~ s(temp, k = 3), 
              data = d,
              spatial = "off",
              spatiotemporal = "off"#,
              #family = student(df = 6)
              )
  
  res <- residuals(m)
  qqnorm(res); qqline(res)
  
  nd <- data.frame(area = a,
                   temp = seq(min(d$temp), max(d$temp), length.out = 50))
  
  nd$pred <- predict(m, newdata = nd)$est
    
  print(ggplot(nd, aes(temp, pred)) +
    geom_point(data = d, aes(temp, k_median), inherit.aes = FALSE) + 
    geom_line())

  preds <- bind_rows(preds, nd)
    
}
#> filter: removed 319 rows (84%), 63 rows remaining

#> filter: removed 328 rows (86%), 54 rows remaining

#> filter: removed 341 rows (89%), 41 rows remaining

#> filter: removed 347 rows (91%), 35 rows remaining

#> filter: removed 334 rows (87%), 48 rows remaining

#> filter: removed 344 rows (90%), 38 rows remaining

#> filter: removed 363 rows (95%), 19 rows remaining

#> filter: removed 348 rows (91%), 34 rows remaining

#> filter: removed 364 rows (95%), 18 rows remaining

#> filter: removed 367 rows (96%), 15 rows remaining

#> filter: removed 370 rows (97%), 12 rows remaining

#> filter: removed 377 rows (99%), 5 rows remaining

Code

  m_all <- sdmTMB(k_mean ~ s(temp, k = 5), 
                  data = VBG,
                  spatial = "off",
                  spatiotemporal = "off"#,
                  #family = student(df = 6)
                  )
  
  res_all <- residuals(m_all)
  qqnorm(res_all); qqline(res_all)

Code

  
  nd_all <- data.frame(temp = seq(min(VBG$temp), max(VBG$temp), length.out = 50))
  
  nd_all$pred <- predict(m_all, newdata = nd_all)$est
    
  ggplot(nd_all, aes(temp, pred)) +
    geom_point(data = VBG, aes(temp, k_median), inherit.aes = FALSE) + 
    geom_line()

Code

preds |>
  ggplot(aes(temp, pred, color = factor(area, order$area))) + 
  geom_point(data = VBG, aes(temp, k_median, color = factor(area, order$area)), size = 0.6) +
  geom_line(aes(temp, pred, group = factor(area)), linewidth = 1) +
  geom_line(data = nd_all, aes(temp, pred),
            linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  #scale_color_manual(values = colors, name = "Area") +
  scale_color_viridis(name = "Area", direction = -1, discrete = TRUE) +
  #facet_wrap(~area) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (C)", y = "Median von Bertalanffy growth coefficient, k") +
  NULL